A soluble model of evolution and extinction dynamics in a rugged fitness landscape. 
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. . . We consider a continuum version of a previously introduced and numerically studied model of 

macroevolution ||^, in which agents evolve by an optimization process in a rugged fitness landscape 
and die due to their competitive interactions. We first formulate dynamical equations for the 
I fitness distribution and the survival probability. Secondly we analytically derive the t~^ law which 

■ characterizes the life time distribution of biological genera. Thirdly we discuss other dynamical 

properties of the model as the rate of extinction and conclude with a brief discussion. 

^ ! 87.10.+e,02.50.f2,05.40.+j,03.20.+i 

T-H ' Aspects of evolution and extinction can be described as emergent behavior in a large set of interacting agents 
-Q, moving stochastically in a rugged fitness landscape |Q. The behavior of the models of Refs. (|-^ stems from 
fluctuations in a time homogeneous stochastic process. This agrees with a commonly held perception, e.g. implied 
when a birth-death process with constant rates Q is used to fit survivorship data and when the size of extinction events 
O i' is presented as a 'kill curve' 0. A quite different paradigm is also frequently met in the literature: Raup and Sepkoski 
Q noted that the apparent decrease of the extinction rate through geological times could be predictable from 
first principles if one argues that general optimization of fitness through evolutionary time should lead to prolonged 
O ^ survival '. Gould |^ uses an unexpected source of statistical data to illustrate evolutionary non-homogeneity as it 
. reveals itself in the 'unreversed, but constantly slowing, improvement in mean fielding average through the history of 
O ' baseball'. Concurring observations from experimental studies of bacterial evolution in a constant environment can be 
c/) found in Ref. |l0|] as well as from numerical experiments on the 'long jump dynamics' of the NK model in Ref. p| . 

In this Letter we consider stochastic evolution in a rugged fitness landscape. The assumptions are the same in 
spirit as those of a previously introduced and numerically studied 'reset' model However, they are here expressed 
in a further simplified way, allowing a (mainly) analytical rather than (mainly) numerical treatment, and leading 
to close form expressions for the survivorship curves and life span distributions, which are of general interest in the 
study of complex evolving systems, biological or not. We use two - somewhat extreme - assumptions in line with a 
^ ' non-stationary evolution paradigm: Firstly, the progeny of individual mutants less fit than the currently dominating 
genotype never establish itself within the population. Then, as a macroscopic evolutionary step can only be triggered 
by a fitness record within the population, the current typical genotype always codes the best solution found 'so far', 
ly-^ ] Secondly, competitive interactions among species depend on fitness in a non-symmetric way, as evolving species only 
• kill their less fit neighbors. The predictions of the present model resemble the behavior of the reset model and are in 
t — ' good agreement with empirical data describing biological genera [|6PJl^-[l5| . 

0^ In the sequel we first derive equations for the fitness distribution of the system P{x,t) and for the probability Wt{T) 

' that a tagged species born at time t survive time r. We then analytically find the dependence of W and the 
O , ensuing dependence of the life-time distribution Rt{T). Next we discuss the parametric t dependence which is 
c/3 ■ not in general analytically available, the effect of averaging over t, and the long time asymptotic behavior of P{x, t) 
J>~>| for different parameter values. We conclude by with a brief assessment of the robustness of the model. 

' To construct a dynamical equation for P{x, t) we proceed in two steps, starting with the limiting case where no 
TT^\ extinctions take place and where, as a consequence of hill climbing in a random fitness landscape, a suitably defined 
t> [p], p^p^ average fitness grows logarithmically: 

><■ D{t)=\og{t + l). (1) 

d ' With no interactions, an initial fitness distribution would be rigidly shifted in (log) time. As D solves the equation 
of motion v{x) = dx/dt — exp(— a;) with initial condition D{Q) = 0, the time evolution of a distribution of non 
interacting agents P{x,t) solves the transport equation: dP{x,t)/dt + d{v{x)P{x,t))/dx = 0. Interactions enter via 
an additional term —gP{x,t)K{P{x,t)), where K is an effective killing rate and where the constant g describes what 
fraction of the system is affected by an evolutionary event. 

Species going extinct vacate a niche, which is refilled at a later time. This in and outflow is expediently accounted 
by introducing a 'limbo' state, which absorbs extinct species, and from which new species emerge at the low fitness 
boundary of the system. A finite upper bound to the total number of species which can coexist implies a conservation 
law : N{t) + P{y, t)dy — 1. With the chosen normalization N{t) is the fraction of species in the limbo state, while 
P{x,t) is the probability density of finding a living species with fitness x. The above considerations lead us to the 
differential equations: 
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-bN{t)+g / P{z,t)K{P)dz (2) 
gP{x,t)K{P) (3) 



dt dx 

where b is the rate at which species are generated at the low fitness end of the system. The corresponding initial and 
boundary conditions are: N{t = 0) = A^o, Vcc : P{x,t = 0) = Po{x), Wt : P{x,t)dx 1 — N{t) < oo, and finally 
yt : P{x = 0,t) = hN{t). 

We consider below a form of the killing rate K which is as close as possible to the reset model: The killing at fitness 
X is taken to depend on the rate of evolutionary change of agents with fitness larger than x: low-fitness agents suffer 
if high fitness agents evolve - but not vice versa. This leads to 

K{v{x)P{x,t)) = {- d{vP)/dx)°'dx = {v{x)P{x,t))", (4) 

J X 

simply expressing the killing rate as the evolutionary current raised to a power. The exponent a just introduced 
allows more generality without unduly complicating the analysis: It accounts in a simplified way for possible (spatial) 
correlations effects in a model where information about individual species is retained. If a < 1 (> 1), a move by 
an old, slowly evolving species triggers a larger (smaller ) cascade of extinctions than one by a young, fast evolving 
species. Figure 1 shows six snaphots of the fitness distribution resulting from the above equations, at times equally 
spaced on a logarithmic scale and for a = 1, 6 = 1 and g — 40. 

A quantity often used to characterize paleontological data is the survivorship curve of a cohort or the closely related 
life span distribution In our treatment the former quantity corresponds to the probability W^t(T) that an agent 
appearing at time t survive time r, while the latter can be found from Wt(r) by differentiation: 

«.M = (5) 

As an agent born at i and alive at time i + t invariably has fitness P>(t) = ln(r + 1) and as the probability of being 
killed in the interval dr is K{P{D{t), t + T))dT, W must obey the differential equation: 

'-}^^-gKiPiDirU + r))r>0, (6) 
dr 

with initial condition Wt(r = 0) = 1. 

Finally, the model extinction rate is simply the fraction of species which die per unit of time, at time t: 

POO 

r{t) =g P{x, t)K{P{x, t))dx = dN/dt + bN. (7) 
Jo 

Note that if 6 — s- oo, then extinct species are immediately replaced, as in Ref. 0]. Furthermore for any b and large bt 
dN/dt is negligible and bN{t) — *■ r{t) so that the extinction closely balances the inflow. 

As a first step towards the solution of Eq.^, we set q = vP and notice that q can be written as q{z{x^t)) with 

dq/dz = ~gq"+\ (8) 

and where z{x,t) satisfies 

dz/dt + v{x)dzldx = l. (9) 

The solution of Eq.^ is simply q — (agz)^^ I . To solve Eq.|| we let A and B be any two functions of a single real 
variable x, which are continuous for x > and which vanish identically for a: < 0. For v = exp(— x), the general solution 
has the form z(a;, t) = e exp(a;) + {l — e)t + A{t+l — exp(a;)) + i?(exp(a;) — (t+l)) for some constant e < 1. Utilizing the 
initial and boundary conditions, we find A{y) — {bN{y))~°', y > and B{y) = {y + l)'^Po{log{y + l))^"—gay, y > 0, 
leading to 

P{x,t)^ T7- (10) 

[gat+ie--t)'^p-"{logie--t))y^" 

for X > D{t), while for x < D{t) we have 
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P{x, t) = — (11) 

[ga{e^ - 1) + {bN{t + 1 - e^))""]'/" 

Note that P is continuous in x, although its derivative will in general be discontinuous at x = D{i)- 

The survival probability of a species born at time t (the survivorship curve of a cohort j^]), can be obtained 
analytically by solving EqJ^. This is so because when inserting D{t) in lieau of x in Eq.^, the t dependence in the 
argument of the (unknown) function N drops out. The solution is: 



Wt{T) 



ibNit))-'^ "1^/" 



_{bN{t))-" + gar 



(12) 



As W vanishes for large r, all species eventually die, regardless of the value of a. This behavior is very desirable from 
a modeling point of view, as it agrees with the fact that by far the largest number of species which ever lived are now 
extinct |14|. The distribution of life spans can be obtained from Eq.|l2| by differentiation, as expressed in Eq.||. If a 



is close to unity, we find a t^^ behavior for Wt(T), and hence a for Rtir), independently of t. 

Averaging these distributions with respect to t over a time window T is needed if the time of appearance of species 
is not precisely known, or if data are scarce. Weighing Rt{T) by the normalized rate at which new species flow into 
the system we obtain: 

Of course, averaging does not change the behavior significantly if T is short compared to the typical lifetime of the 
species. In the opposite limit, the behavior is also maintained if N{t) does not not vanish 'too rapidly' in the limit 
t ^ oo. To better appreciate this last point, we use Eq.|l^ in conjunction with Eqsjl^ and ^, and express A^~" by 
Wt{T)^ obtaining: 

, ,1 fJ^""(l-M^°(T))l + l/"dy 

Jo bN{y)dy 

Even though this integral cannot be evaluated explicitly, Eq.l^ shows that the r dependence of the integrand is 
negligible if the inequality agr > {bN{y))~°' holds throughout the integration interval. The r dependence stemming 
from the limits of the integrals can also be ignored for t << T. Hence R oc ^--i-i/"^ similarly to the non-averaged 
case. As shown later, when a close to unity and g sufficiently large, the model yields r{t) w bN{t) > , with 5 
close to 0.5, which means that even though t << T the relation r > T^" can be fulfilled. 

We now restrict ourselves to a limiting case in which N(t = 0) = 1 which is formally at variance with our boundary 
conditions. However, a limiting process shows that the relevant expression for P{x,t) for x < D{t) remains EqHlL 
while P — Q ior X > D{t). A non- linear equation for N{t) is now obtained by integration of Eq.^ followed by a 
change of variables. The result is 

1 - N{t) = f ^ (15) 

io [ga[t-y) + {bN{y)r<-f''' 

Differentiating Eq.^ with respect to time, and utilizing Eq.^, we find the extinction rate: 

= /* '-^ ^ (16) 

Jo [5a(i-y) + (6iV(y))-"]^+^/" 

A closed form solution of these (equivalent) integral equations could not be found in the general case. We notice 
however a major difference in the asymptotic behavior for a < 1 and a > 1. In both cases the time independent 
function Poo{x) obtained by taking i — > cx) and by setting N{t) = Noo ^ in Eq.]!!] formally satisfies the model 
equations. However, only for a < 1 is Poo{x) normalizable and thus a true solution. The corresponding steady state 
value of N , Noo is then implicitly given by the relation 1 — Noo = {bNoo)^~°' / {g{^ — a)), which always has a solution 
in the unit interval. 

In the case a > 1, normalizability of P{x,t) requires that N{t) ^ for t — > c». No steady state solution can then 
exists, since P{x, t) vanishes with t at any fixed a:, as e.g. in the familiar case of simple diffusion on the infinite line. 
For a < 1, the steady state solution is strictly speaking only approached logarithmically due to the form of D{t). 
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Neglecting this logarithmic corrections we see from Fig. 2 a power-law approach to a quasistationary behavior over a 
substantial time range. 

For long times a = 1 and bt >> 1 the term dN/dt in Eq.||is negligible, and bN{t) w r(t). In this limit we can also 
neglect N compared to one, thus finding the following approximate equation for r{t): 1 = g(^t-y)+r{y)-^ which has 
the solution r{t) — (gt)^^. 

Fig. 2 shows a r{t) vs. t for a — 0.95 b — 1 and several g values. As noted, for a wide time span, r cx 
where 7 decreases with increasing g, similarly to the result obtained in the the simulations of the reset model [Q. No 
qualitative changes are observed when varying a in a small range below one, or when changing b. In summary, for 
a slightly below one, and g is sufficiently large, the life span distribution (averaged or not) decays algebraically with 
an exponent slightly above —2 and the rate of extinction decays with an exponent close to —0.5 untill it reaches a 
regime of hardly detectable change. 

The most comprehensive empirical life span distribution available, comprising about 17500 extinct genera of marine 
animals has been tabulated by Raup |^] from data compiled by Sepkoski These data cover about 100 million 

years and display a very clear dependence in a log-log plot [Qjsj over this range, which concurs with the behavior 
of our Rt{T). More recent analysis by Baumiller |]l5| of several data set describing crinoid survivorship - our Wt{T) - 
over a comparable time span in part concur with a t~^ law, and hence with a t^^ law for the life-time distribution. 
Finally, survivorship curves for european mammals were considered by Stanley [^2|. These data span approximately 
3 million years stretching to the Wiirm period and include much fewer species. The distribution of lifetimes deviates 
from a t~'^ law by having an extra 'hump' approximately in the middle of the time range. 

Paleontological data are commonly interpreted using a birth and death model [^|5|, in which non- interacting 
species are born and die with two distinct constant rates of speciation and extinction, A and /i, and where the genus 
becomes extinct together with its last species. Interestingly, the survivorship formula generated by this model is, for 
X = fj, and for an initial number of species equal to one, identical to our Eq.l^ - with a = 1, as far as its dependence 
on the life time -our r- goes. By continuity so arc the model predictions in the often recurring situation when X ^ fi. 
The similarity in the formulae is however contingent to the initial condition and should be regarded as accidental [ p^ . 

In line with the conclusion of the reset model, we have shown analytically that a large body of data describing 
evolution on coarse scales of time and taxonomical level can be explained by two very simple ideas: 1) that fitness 
records in random searching trigger evolutionary events, and 2) that the species competition is 'asymmetric', with 
high fitness species being more resilient. 

The robustness of this approach has already been analyzed to some extent: The effect of additional and extern ally 
imposed random killings of a fraction of the agents ~ mimicking catastrophies - has been studied by M. Brandt pi[ , 
who found that the life-span distribution was not affected. This is to be expected, as even very large mass extinction 
events - in the model as well as in reality - only account for a small fraction of all extinctions. We also explored 
other choices of the killing term in Eq.^, finding that the l/t^ law disappears if the asymmetry of the interspecies 
interactions is removed, with the possible exception of special values of the coupling constants. 

After this paper was submitted the author became aware of a preprint by Manrubia and Paczuski , which also 
treats evolution and extinctions by means of a transport equation, an finds a t^"^ life-time distribution, albeit the 
basic dynamical mechanism is quite different from ours. 
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FIG. 1. For 1 < n < 6 the n'th plot in the figure depicts the fitness distribution P{x,t) at time t = 10" The parameters 
used are a = 1, 6 = 1 and g = 40. 




FIG. 2. The numerically obtained extinction rate is plotted vs time, for b = 1, a = 0.95 and for g = 5 (o), 10 (.), 20 (*), 
30 (+) and 40 (x). The decay is approximately a power-law with exponents equal to —0.89, —0.85, —0.75, —0.69 and —0.56. 
These slopes are visualized by the full lines. 
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